package com.clei.utils;

import com.clei.dto.geo.BoundingBox;
import com.clei.dto.geo.Point;

/**
 * 坐标处理工具类
 * <p>
 * WGS84坐标系：即地球坐标系，国际上通用的坐标系。设备一般包含GPS芯片或者北斗芯片获取的经纬度为WGS84地理坐标系。
 * 谷歌地图采用的是WGS84地理坐标系（中国范围除外,谷歌中国地图采用的是GCJ02地理坐标系。)
 * <p>
 * GCJ02坐标系：即火星坐标系，WGS84坐标系经加密后的坐标系。
 * 出于国家安全考虑，国内所有导航电子地图必须使用国家测绘局制定的加密坐标系统，即将一个真实的经纬度坐标加密成一个不正确的经纬度坐标。
 * <p>
 * BD09坐标系：即百度坐标系，GCJ02坐标系经加密后的坐标系。搜狗坐标系、图吧坐标系等，估计也是在GCJ02基础上加密而成的。
 * <p>
 * 高德MapABC地图API  火星坐标
 * 腾讯搜搜地图API  火星坐标
 * 阿里云地图API  火星坐标
 * 灵图51ditu地图API  火星坐标
 * <p>
 * 百度地图API  百度坐标
 * 搜狐搜狗地图API  搜狗坐标
 * 图吧MapBar地图API  图吧坐标
 *
 * @author Y
 * @date 2022-05-07
 */
public final class CoordinateUtil {


    /**
     * 圆周率π
     */
    private static final double PI = 3.1415926535897932384626D;

    /**
     * 火星坐标系与百度坐标系转换的中间量
     */
    private static final double X_PI = 3.14159265358979324 * 3000.0 / 180.0D;

    // Krasovsky 1940
    // 长半轴a = 6378245.0, 1/f = 298.3
    // b = a * (1 - f)
    // 扁率ee = (a^2 - b^2) / a^2;

    /**
     * 长半轴，克拉索夫斯基椭球参数
     */
    private static final double SEMI_MAJOR = 6378245.0D;

    /**
     * 扁率，克拉索夫斯基椭球参数第一偏心率平方
     */
    private static final double FLATTENING = 0.00669342162296594323D;

    /**
     * 赤道半径 米
     */
    private static final double EQUATOR_RADIUS = 6378137;

    /**
     * 地球平均半径 米
     */
    private static final double EARTH_AVG_RADIUS = 6371393;

    /**
     * WGS84=>GCJ02 地球坐标系=>火星坐标系
     *
     * @param point 地球坐标
     * @return 火星坐标
     */
    public static Point wgs84ToGcj02(Point point) {
        return wgs84ToGcj02(point.getLongitude(), point.getLatitude());
    }

    /**
     * WGS84=>GCJ02 地球坐标系=>火星坐标系
     *
     * @param longitude 经度
     * @param latitude  纬度
     * @return 火星坐标
     */
    public static Point wgs84ToGcj02(double longitude, double latitude) {
        if (outOfChina(longitude, latitude)) {
            return new Point(longitude, latitude);
        }

        double[] offset = offset(longitude, latitude);
        double mgLongitude = longitude + offset[0];
        double mgLatitude = latitude + offset[1];

        return new Point(mgLongitude, mgLatitude);
    }

    /**
     * GCJ02=>WGS84 火星坐标系=>地球坐标系(粗略)
     *
     * @param point 火星坐标
     * @return 地球坐标
     */
    public static Point gcj02ToWgs84(Point point) {
        return gcj02ToWgs84(point.getLongitude(), point.getLatitude());
    }

    /**
     * GCJ02=>WGS84 火星坐标系=>地球坐标系(粗略)
     *
     * @param longitude 经度
     * @param latitude  纬度
     * @return 地球坐标
     */
    public static Point gcj02ToWgs84(double longitude, double latitude) {
        if (outOfChina(longitude, latitude)) {
            return new Point(longitude, latitude);
        }

        double[] offset = offset(longitude, latitude);
        double mgLongitude = longitude - offset[0];
        double mgLatitude = latitude - offset[1];

        return new Point(mgLongitude, mgLatitude);
    }

    /**
     * GCJ02=>WGS84 火星坐标系=>地球坐标系（精确）
     *
     * @param point 火星坐标
     * @return 地球坐标
     */
    public static Point gcj02ToWgs84Exactly(Point point) {
        return gcj02ToWgs84Exactly(point.getLongitude(), point.getLatitude());
    }

    /**
     * GCJ02=>WGS84 火星坐标系=>地球坐标系（精确）
     *
     * @param longitude 经度
     * @param latitude  纬度
     * @return 地球坐标
     */
    public static Point gcj02ToWgs84Exactly(double longitude, double latitude) {
        if (outOfChina(longitude, latitude)) {
            return new Point(longitude, latitude);
        }

        double initDelta = 0.01;
        double threshold = 0.000000001;
        double dLatitude = initDelta, dLongitude = initDelta;
        double mLatitude = latitude - dLatitude, mLongitude = longitude - dLongitude;
        double pLatitude = latitude + dLatitude, pLongitude = longitude + dLongitude;
        double wgsLatitude, wgsLongitude, i = 0;
        while (true) {
            wgsLatitude = (mLatitude + pLatitude) / 2;
            wgsLongitude = (mLongitude + pLongitude) / 2;
            Point point = wgs84ToGcj02(wgsLongitude, wgsLatitude);
            dLongitude = point.getLongitude() - longitude;
            dLatitude = point.getLatitude() - latitude;
            if ((Math.abs(dLatitude) < threshold) && (Math.abs(dLongitude) < threshold)) {
                break;
            }

            if (dLatitude > 0) {
                pLatitude = wgsLatitude;
            } else {
                mLatitude = wgsLatitude;
            }
            if (dLongitude > 0) {
                pLongitude = wgsLongitude;
            } else {
                mLongitude = wgsLongitude;
            }

            if (++i > 10000) {
                break;
            }
        }

        return new Point(wgsLongitude, wgsLatitude);
    }

    /**
     * GCJ-02=>BD09 火星坐标系=>百度坐标系
     *
     * @param point 火星坐标
     * @return 百度坐标
     */
    public static Point gcj02ToBd09(Point point) {
        return gcj02ToBd09(point.getLongitude(), point.getLatitude());
    }

    /**
     * GCJ-02=>BD09 火星坐标系=>百度坐标系
     *
     * @param longitude 经度
     * @param latitude  纬度
     * @return 百度坐标
     */
    public static Point gcj02ToBd09(double longitude, double latitude) {
        double z = Math.sqrt(longitude * longitude + latitude * latitude) + 0.00002 * Math.sin(latitude * X_PI);
        double theta = Math.atan2(latitude, longitude) + 0.000003 * Math.cos(longitude * X_PI);
        double bdLongitude = z * Math.cos(theta) + 0.0065;
        double bdLatitude = z * Math.sin(theta) + 0.006;
        return new Point(bdLongitude, bdLatitude);
    }

    /**
     * BD09=>GCJ-02 百度坐标系=>火星坐标系
     *
     * @param point 百度坐标
     * @return 火星坐标
     */
    public static Point bd09ToGcj02(Point point) {
        return bd09ToGcj02(point.getLongitude(), point.getLatitude());
    }

    /**
     * BD09=>GCJ-02 百度坐标系=>火星坐标系
     *
     * @param longitude 经度
     * @param latitude  纬度
     * @return 火星坐标
     */
    public static Point bd09ToGcj02(double longitude, double latitude) {
        double x = longitude - 0.0065;
        double y = latitude - 0.006;
        double z = Math.sqrt(x * x + y * y) - 0.00002 * Math.sin(y * X_PI);
        double theta = Math.atan2(y, x) - 0.000003 * Math.cos(x * X_PI);
        double gcjLongitude = z * Math.cos(theta);
        double gcjLatitude = z * Math.sin(theta);
        return new Point(gcjLongitude, gcjLatitude);
    }

    /**
     * WGS84=>BD09 地球坐标系=>百度坐标系
     *
     * @param point 地球坐标
     * @return 百度坐标
     */
    public static Point wgs84ToBd09(Point point) {
        return wgs84ToBd09(point.getLongitude(), point.getLatitude());
    }

    /**
     * WGS84=>BD09 地球坐标系=>百度坐标系
     *
     * @param longitude 经度
     * @param latitude  纬度
     * @return 百度坐标
     */
    public static Point wgs84ToBd09(double longitude, double latitude) {
        Point point = wgs84ToGcj02(longitude, latitude);
        return gcj02ToBd09(point.getLongitude(), point.getLatitude());
    }

    /**
     * BD09=>WGS84 百度坐标系=>地球坐标系
     *
     * @param point 百度坐标
     * @return 地球坐标
     */
    public static Point bd09ToWgs84(Point point) {
        return bd09ToWgs84(point.getLongitude(), point.getLatitude());
    }

    /**
     * BD09=>WGS84 百度坐标系=>地球坐标系
     *
     * @param longitude 经度
     * @param latitude  纬度
     * @return 地球坐标
     */
    public static Point bd09ToWgs84(double longitude, double latitude) {
        Point point = bd09ToGcj02(longitude, latitude);
        return gcj02ToWgs84(point.getLongitude(), point.getLatitude());
    }

    /**
     * 判断是否在中国境外
     *
     * @param point 地球坐标
     * @return 中国境外返回true, 境内返回false
     */
    public static boolean outOfChina(Point point) {
        return outOfChina(point.getLongitude(), point.getLatitude());
    }

    /**
     * 判断是否在中国境外
     *
     * @param longitude 经度
     * @param latitude  纬度
     * @return 中国境外返回true, 境内返回false
     */
    public static boolean outOfChina(double longitude, double latitude) {
        if (longitude < 72.004 || longitude > 137.8347) {
            return true;
        }
        if (latitude < 0.8293 || latitude > 55.8271) {
            return true;
        }
        return false;
    }

    /**
     * 坐标系转换
     * EPSG4326(WGS84) -> EPSG3857
     *
     * @param longitude 经度
     * @param latitude  纬度
     * @return EPSG3857坐标点
     */
    public static double[] epsg4326To3857(double longitude, double latitude) {
        double earthRad = 6378137.0;
        double x = longitude * Math.PI / 180 * earthRad;
        double a = latitude * Math.PI / 180;
        double y = earthRad / 2 * Math.log((1.0 + Math.sin(a)) / (1.0 - Math.sin(a)));
        return new double[]{x, y};
    }

    /**
     * 两个经纬度之间的距离 米
     *
     * @param lon1 经度1
     * @param lat1 纬度1
     * @param lon2 经度2
     * @param lat2 纬度2
     * @return 距离 米
     */
    public static double distance(double lon1, double lat1, double lon2, double lat2) {
        double radLon1 = Math.toRadians(lon1);
        double radLat1 = Math.toRadians(lat1);
        double radLon2 = Math.toRadians(lon2);
        double radLat2 = Math.toRadians(lat2);

        double latDiff = radLat1 - radLat2;
        double lonDiff = radLon1 - radLon2;
        double distance = 2 * Math.asin(Math.sqrt(Math.pow(Math.sin(latDiff / 2), 2)
                + Math.cos(radLat1) * Math.cos(radLat2)
                * Math.pow(Math.sin(lonDiff / 2), 2)));
        distance = distance * EQUATOR_RADIUS;
        return distance;
    }

    /**
     * 两个经纬度点之间的距离 米
     *
     * @param p1 点1
     * @param p2 点2
     * @return 距离 米
     */
    public static double distance(Point p1, Point p2) {
        return distance(p1.getLongitude(), p1.getLatitude(), p2.getLongitude(), p2.getLatitude());
    }

    /**
     * 两个经纬度之间的距离 米
     * 高德地图计算方法
     *
     * @param lon1 经度1
     * @param lat1 纬度1
     * @param lon2 经度2
     * @param lat2 纬度2
     * @return 距离 米
     */
    public static double gdDistance(double lon1, double lat1, double lon2, double lat2) {
        if (lon1 == 0 || lat1 == 0 || lat2 == 0 || lon2 == 0) {
            return -1.0;
        }
        lon1 *= 0.01745329251994329;
        lat1 *= 0.01745329251994329;
        lon2 *= 0.01745329251994329;
        lat2 *= 0.01745329251994329;
        double var1 = Math.sin(lon1);
        double var2 = Math.sin(lat1);
        double var3 = Math.cos(lon1);
        double var4 = Math.cos(lat1);
        double var5 = Math.sin(lon2);
        double var6 = Math.sin(lat2);
        double var7 = Math.cos(lon2);
        double var8 = Math.cos(lat2);
        double[] var10 = new double[3];
        double[] var20 = new double[3];
        var10[0] = var4 * var3;
        var10[1] = var4 * var1;
        var10[2] = var2;
        var20[0] = var8 * var7;
        var20[1] = var8 * var5;
        var20[2] = var6;

        return Math.asin(Math.sqrt((var10[0] - var20[0]) * (var10[0] - var20[0]) + (var10[1] - var20[1]) * (var10[1] - var20[1]) + (var10[2] - var20[2]) * (var10[2] - var20[2])) / 2.0) * 1.27420015798544E7;
    }

    /**
     * 两个经纬度之间的距离 米
     * 高德地图计算方法
     *
     * @param p1 点1
     * @param p2 点2
     * @return 距离 米
     */
    public static double gdDistance(Point p1, Point p2) {
        return gdDistance(p1.getLongitude(), p1.getLatitude(), p2.getLongitude(), p2.getLatitude());
    }

    /**
     * 以点位为中心 边长为2 * distance的正方形边框
     *
     * @param lon      经度
     * @param lat      纬度
     * @param distance 距离 米
     * @return BoundingBox
     */
    public static BoundingBox boundingBox(double lon, double lat, double distance) {
        double dLon = 2 * Math.asin(Math.sin(distance / (2 * EARTH_AVG_RADIUS)) / Math.cos(Math.toRadians(lat)));
        double dLat = distance / EARTH_AVG_RADIUS;
        double minLon = lon - Math.toDegrees(dLon);
        double maxLon = lon + Math.toDegrees(dLon);
        double minLat = lat - Math.toDegrees(dLat);
        double maxLat = lat + Math.toDegrees(dLat);
        return new BoundingBox(minLon, minLat, maxLon, maxLat);
    }


    /**
     * 经度偏移量
     *
     * @param longitude 经度
     * @param latitude  纬度
     * @return 转换后的经度
     */
    private static double transformLongitude(double longitude, double latitude) {
        double ret = 300.0 + longitude + 2.0 * latitude + 0.1 * longitude * longitude + 0.1 * longitude * latitude + 0.1 * Math.sqrt(Math.abs(longitude));
        ret += (20.0 * Math.sin(6.0 * longitude * PI) + 20.0 * Math.sin(2.0 * longitude * PI)) * 2.0 / 3.0;
        ret += (20.0 * Math.sin(longitude * PI) + 40.0 * Math.sin(longitude / 3.0 * PI)) * 2.0 / 3.0;
        ret += (150.0 * Math.sin(longitude / 12.0 * PI) + 300.0 * Math.sin(longitude / 30.0 * PI)) * 2.0 / 3.0;
        return ret;
    }

    /**
     * 纬度偏移量
     *
     * @param longitude 经度
     * @param latitude  纬度
     * @return 转换后的纬度
     */
    private static double transformLatitude(double longitude, double latitude) {
        double ret = -100.0 + 2.0 * longitude + 3.0 * latitude + 0.2 * latitude * latitude + 0.1 * longitude * latitude + 0.2 * Math.sqrt(Math.abs(longitude));
        ret += (20.0 * Math.sin(6.0 * longitude * PI) + 20.0 * Math.sin(2.0 * longitude * PI)) * 2.0 / 3.0;
        ret += (20.0 * Math.sin(latitude * PI) + 40.0 * Math.sin(latitude / 3.0 * PI)) * 2.0 / 3.0;
        ret += (160.0 * Math.sin(latitude / 12.0 * PI) + 320 * Math.sin(latitude * PI / 30.0)) * 2.0 / 3.0;
        return ret;
    }

    /**
     * 偏移量
     *
     * @param longitude 经度
     * @param latitude  纬度
     * @return 经纬度偏移量
     */
    private static double[] offset(double longitude, double latitude) {
        double[] lngLat = new double[2];
        double dLongitude = transformLongitude(longitude - 105.0, latitude - 35.0);
        double dLatitude = transformLatitude(longitude - 105.0, latitude - 35.0);
        double radLatitude = latitude / 180.0 * PI;
        double magic = Math.sin(radLatitude);
        magic = 1 - FLATTENING * magic * magic;
        double sqrtMagic = Math.sqrt(magic);
        dLongitude = (dLongitude * 180.0) / (SEMI_MAJOR / sqrtMagic * Math.cos(radLatitude) * PI);
        dLatitude = (dLatitude * 180.0) / ((SEMI_MAJOR * (1 - FLATTENING)) / (magic * sqrtMagic) * PI);
        lngLat[0] = dLongitude;
        lngLat[1] = dLatitude;
        return lngLat;
    }
}